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We review the problem of finding an apparent horizon in Cauchy data (E, g a b, K a b) in three space 
dimensions without symmetries. We describe a family of algorithms which includes the pseudo- 
spectral apparent horizon finder of Nakamura et al. and the curvature flow method proposed by 
Tod as special cases. We suggest that other algorithms in the family may combine the speed of the 
former with the robustness of the latter. A numerical implementation for Cauchy data given on a 
grid in Cartesian coordinates is described, and tested on Brill-Lindquist and Kerr initial data. The 
new algorithm appears faster and more robust than previous ones. 



I. INTRODUCTION 



An important task in numerical relativity is locating black holes in numerically generated spacetimes, both for 
technical purposes and for extracting physical information. A black hole is a region of spacetime out of which no null 
geodesies escape to infinity. The boundary of the black hole, the event horizon, is formed by those outward-going, 
' future-directed null geodesies which neither fall into the singularity nor escape to null infinity. The event horizon 
contains important geometric information about the spacetime. It is a global construction and can in principle only 
be determined when the entire spacetime is known. In practice, one can obtain a good approximation to the event 
I/"") \ horizon within a finite spacetime region, once the black hole has settled down to a stationary state. By definition, 
the event horizon repels future-directed null geodesies, but attracts past ones. One can then evolve past-directed null 
geodesies back through the spacetime, and find the event horizon as the surface to which they are attracted [jij. 

Locating black holes is crucial in numerical relativity also for a technical reason: Spacetime slicings which avoid 
black holes rapidly become singular. Instead one would like to excise a spacetime region just inside the event horizon 
from the numerical domain during the numerical evolution, using the fact that it cannot influence events outside the 
' black hole. During the time evolution, however, one does not yet know where the event horizon is. Instead one needs 
I | to use the poor man's event horizon, the apparent horizon. 

J-Hi An apparent horizon (AH) is defined within a single time slice, or spacelike hypersurface S, namely as a smooth 
embedded 2-surface whose outgoing normal null geodesies have zero expansion. There may be one such surface 
^ . enclosing another one, in which case the outermost one is the apparent horizon. If one combines the apparent horizon 
on each time slice into a 3-dimensional surface, this world tube will depend on the slicing, and can be discontinuous. 
Nevertheless one can show that if an apparent horizon exists on a given time slice, it must be inside a black hole ||. 
The converse is not true: there are slicings of black hole spacetimes without any apparent horizons ||. For numerical 
purposes one simply hopes that this case is unusual and that the apparent horizon gives a reasonable indication of 
the event horizon. 

A wide variety of numerical algorithms for finding AHs have been explored or suggested. For the purpose of excising 
the black hole region, one needs to find the apparent horizon frequently, perhaps at each time step. When two black 
holes collide, a new AH enveloping the two separate ones appears suddenly. Therefore the main requirements are 
speed and finding the AH from scratch, without a good initial guess. Precision is less important for black hole excision, 
although a safe error estimate is, so that one can be sure not to excise too much and inject unphysical boundary 
conditions. 

In spherical symmetry, the AH problem reduces to an algebraic equation. In axisymmetry, it reduces to an ordinary 
differential equation with periodic boundary conditions. In this paper we shall be concerned exclusively with the 3- 
dimensional problem without any symmetries, either continuous or discrete, where one deals with a highly nonlinear 
elliptic problem on a closed 2-surface. In practice this will always be the 2-sphere, or several disconnected 2-spheres 
H , which can be treated separately. 
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All 3D AH finder algorithms proposed so far can be classified according to a few key choices, which can be made 
independently one from another. How are candidate AHs represented? One can parameterize an embedded 2-surface 
either by introducing coordinates on it, or as a level set of a function on the 3-dimensional slice in which it is embedded. 
How is the curvature of the candidate AH calculated? One can discretize the necessary spatial derivatives by finite 
differencing, finite elements, or pseudo-spectral methods. A third fundamental choice is between solving the elliptic 
problem directly, or converting it into a parabolic problem, in which the solution of the elliptic problem is approached 
during an evolution in an unphysical time parameter. The distinction between these last two approaches is not sharp 
in practice. On the one hand one always solves a nonlinear elliptic problem by iteration. On the other, numerical 
implementation of any parabolic approach requires an implicit "time" step for stability, thus posing a new elliptic 
problem that becomes equivalent to the original one in the limit of an infinitely large time step. We now discuss 
previous AH finders in terms of these choices. 

Nakamura, Koshima and Oohara || represent the AH in spherical coordinates as r — h(6,tp). We note that this 
requires the surface to be a 2-sphere, and star-shaped (convex) around the point r = 0. The shape function h(9, ip) is 
expanded in spherical harmonics. This spectral decomposition is used to calculate the derivatives of h in formulating 
the elliptic problem. The orthonormality and completeness of the spherical harmonics is used to subtract the linear 
elliptic operator L 2 from the nonlinear elliptic problem and invert it. This gives rise to an iteration prescription. We 
shall see that this iteration can also be described as the discretization in unphysical time of a parabolic problem. It 
remains unclear from H by what method the spectral decomposition back and forth is carried out for Cauchy data 
which arc only known in numerical form and on a grid. The Nakamura et al. algorithm has been independently coded 
and tested, and extended in various directions by Kemball and Bishop Q . They report exponential convergence, good 
robustness, and high precision unless the point r = is close to the AH. 

Tod |?J has proposed a geometrically defined flow under which a trial 2-surface evolves to the AH. For time-symmetric 
slices, the AH problem reduces to that of finding a minimal surface, and Tod's prescription to mean curvature flow. 
This is well-known to converge to minimal surfaces. On non time-symmetric slices, only lower order terms are added 
to the problem, so that one may hope that Tod's flow also converges for such data in practice. Tod's algorithm is 
parabolic, without specifying how the surface is represented or differenced. Tod's algorithm has been implemented 
numerically by Bernstein J^] using finite differencing in coordinates introduced on the surface. He discusses stable 
extrinsic algorithms for parabolic problems, and reports good results in axisymmetry using one of them, but technical 
problems to do with finite differencing on the sphere in the general case. Pasch || has implemented mean curvature 
flow representing the test surface as a level set. This allows for a change of topology during the evolution. He has 
successfully tested the algorithm using Brill-Lindquist data for 1, 2 or 3 black holes, using a fast implicit time evolution 
package, and finite differencing on a Cartesian grid in the embedding space. 

Thornburg JTcfl attacks the elliptic problem directly using finite differencing on a square (9, ip) grid, and Newton's 
method to solve the discretized equations. He calculates the Jacobian required for Newton's method by first linearizing 
the differential equations, then finite differencing the result. This is more efficient than numerical differentiation. He 
finds high precision results, but a nonlinear instability in high-frequency modes. Huq has extended Newton's 
method to data without symmetries in Cartesian-type coordinates. 

The NCSA/WashU algorithm jl2],[l3| uses the parameterization r — h(9, ip), and a spectral decomposition to param- 
eterize h and calculate its derivatives. The discretized elliptic problem is solved by applying a standard minimization 
algorithm to the sum of H 2 over all surface points. The spectral basis is not required to be orthonormal for this 
purpose. Baumgarte et al. [[TiJ have implemented the NCSA/WashU algorithm independently, with the difference 
that they use the true spherical harmonics as a basis. Both algorithms locate points on the 2-surface on a square 

ip) grid, interpolating data from the 3-dimensional Cartesian grid used in the 3+1 time evolution. 

In this paper we review different previous algorithms in one common, fully covariant notation. This analysis 
suggests to us a new algorithm which combines essential ideas of the algorithms of Tod and Nakamura et al. From our 
analysis we expect this algorithm to be as fast as that of Nakamura et al. (and therefore much faster than existing 
implementations of Tod's algorithm) , while being as robust in practice as that of Tod. We describe the details of a 
numerical implementation of this algorithm, and some initial tests. The results are encouraging. 

The paper is organized as follows. In section II we set up the mathematical formalism of the problem. We begin 
by deriving the differential equation that determines an apparent horizon in II. A. In II. B we discuss different ways of 
parameterizing apparent horizon candidates, that is, smooth embedded 2-surfaces, and in II. C we provide tools for 
spectral methods on the 2-sphere. In section III we review various algorithms for finding AHs, namely the pseudo- 
spectral algorithm of Nakamura et al. |5), Jacobi's method, and the generalized mean curvature flow suggested by 
Tod @. We then build on this review by presenting a family of algorithm which contains the previous algorithms 
as limiting cases, and suggesting that in the middle of the family there are algorithms that perform better than the 
limiting members. In section IV we describe a numerical implementation of our proposed algorithm. In section V we 
test its performance in finding apparent horizons in Brill-Lindquist and Kerr data given in Cartesian coordinates on 
a grid. 
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II. MATHEMATICAL PRELIMINARIES 



A. The apparent horizon equation 

Here we give a brief derivation (see also of the differential equation that we try to solve in the remainder of 
the paper, both to give a complete presentation of the problem and to fix notation. Throughout the paper, lower-case 
Latin indices from the beginning of the alphabet indicate abstract index notation. Indices from the middle of the 
alphabet indicate 3-dimensional tensor components. Our signature convention is ( — h H — h). 

We begin with a series of definitions. Let (M, ^g a b) be a spacetime, and Vi the covariant derivative associated 
with ( 4 'g a b- (We use this notation to reserve the symbols g a b and V Q for 3 dimensions.) Let E be a smooth spacelike 
hypersurface, and let n a be the future-pointing unit timelike normal to E. Then ^ g a b gives rise to Cauchy data 

gab = {i) gab + n a n b , K ab = ~g a c W{ 4) n b = -V a n b , (1) 

on E, where g ab is the positive definite 3-metric induced on E and K ab is the extrinsic curvature of E. V a is the 
covariant derivative associated with g ab . Let <S* be a closed smooth hypersurface of E, which means it is two-dimensional 
and spacelike, and s a its unit outward pointing normal in E, which is also spacelike, and normal to n a . g ab induces a 
positive definite 2-metric 

m ab = g a b - s a Sb = (4) .g Q fc + n a n b - s a s b (2) 

on S. Let k a be the future-pointing null geodesic congruence whose projection on E is orthogonal to S, that is 

PVi 4 'fc^0, fc a fc a = 0, m ab k a \ s = Q. (3) 

Then k a describes light rays leaving S normally from the point of view of an observer whose instantaneous simultaneity 
is E. Clearly k a depends not only on the spacetime and on S, but also on E. Let H be the expansion of that congruence, 

H = X7^k a . (4) 

We would like to express H in terms of the Cauchy data (E, g ab , K ab ), and the normal s a to S. The crucial step is 
to note that, up to an overall factor, 

k a \ s =s a +n a . (5) 

Clearly this obeys the conditions m ab k a = and k a k a = on S. We continue k a away from S by the remaining 
condition k aT Ua^k b = 0. We also continue s a away from S in E assuming that it retains unit length, but otherwise in 
an arbitrary manner. Then we have, on S, 

H= ( V fc Vl 4 )fc b = (g ab - n a n b )V^k b 

= g ah ^ ( a\sb + n b ) - (k a - s a )(k b - s^V^h 

= g ab ^s b +g ab ^n b - (k b - s b )[k a V^k b ] + s a [k b V^k b ] - s^^s,} - S a s b V^n b 

= V a s a - K + s a s b K ab , (6) 

where K = ^ g ab K ab — g ab K ab is the trace of the extrinsic curvature. All terms in square brackets vanish individually 
by definition. 

A smooth embedded closed surface with outward pointing unit normal s a that obeys H = everywhere on S is 
called a marginally outer trapped surface. The outermost of such surfaces, if one or more exist in E, is called the 
apparent horizon in E On the one hand this definition is global in E, which makes finding an apparent horizon a 
non-trivial problem. On the other, it is local in time, as H depends only on the Cauchy data (g a b, Kab) on a single 
slice E. If one fixes the slicing of a given spacetime, calculates the apparent horizon on each slice, and then combines 
the apparent horizons on each slice to obtain a timelike, 2 + 1 dimensional world-tube, this world-tube depends on the 
slicing. This is in contrast to the event horizon, which depends globally on the entire spacetime, but is independent 
of the slicing. 
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B. Characterizing closed two-surfaces 



Before we can discuss solving the apparent horizon equation H = in practice, we need to parameterize candidate 
apparent horizons, that is, two-dimensional, smooth, closed surfaces S embedded in S. 

Let x l be coordinates on E. One way of parameterizing S is then to introduce coordinates £ A on S (at least locally), 
and give a map x l = X l (£ ). In this case, the topology of S is fixed in advance. Furthermore, different functions X % 
describe the same abstract surface S, corresponding to a change of coordinates £ on S. 

A different way of parameterizing S is as a level set 

Fix') = 0. (7) 

As long as the form of F(x l ) is not restricted, this has the advantage of S allowing to have arbitrary topology. In 
particular, S can be disconnected. Again, many different functions F(x l ) describe the same abstract surface S, as 
long as they have the one level set F(x l ) = in common. 

It is straightforward to express H as a function of F and its derivatives. The unit normal (with respect to the 
3-metric g ab ) of any level set of F is 

s a = |VF|-V 6 V 6 F, where | VF| = [g ab V a FV b F) 1/2 . (8) 

Direct substitution now gives 

H = (g ab - |VFr 2 V FV 6 F) QVFl^VaVftF - K ab ) = m ab (\VF\^V a V b F - K ab ) . (9) 

H is therefore a quasi-linear second order differential operator acting on F. 

Now we come back to the problem that different functions F(x,y,z) describe the same abstract surface S. A 
possible gauge condition would be to make F harmonic with respect to a background metric g ab , or with respect to 
the physical 3-mctric g ab . Then its value everywhere depends only on its value on a suitable two-dimensional surface, 
such as the boundary of the numerical domain. Here, instead, we follow several previous authors in restricting F to 
the form 

F(x l ) = r(x') - h[9(x l ) : p(x% (10) 

where (r, 9, ip) are related to a set of Cartesian coordinates x % in the usual way, namely x = r sin 9 cos ip, y — r sin 9 sin tp 
and z = rcos9. The overall sign of F has been chosen so that s a given in (||) points outward. This parameterization 
is equivalent to X l (9, <p) = x l [r = h(9, ip), 9, tp]. The obvious disadvantages of restricting F to this form are that the 
topology of S must be S 2 , and that S must be star-shaped around the coordinate origin r — 0. The advantages are 
that surfaces S correspond uniquely to functions h, and that we can use the natural basis {ll m } for expanding the 
function h. 

Considered as a quasi- linear differential operator acting on F(x l ), H is not elliptic in three dimensions, because 
one of the three eigenvalues of m a b , the one with eigenvector s a , is zero. Considered as a differential operator in two 
dimensions acting on h(9, ip), it is elliptic. In this two-dimensional interpretation it is nonlinear not only through the 
explicit appearance of W a F in the coefficients of W a ^7bF, but also through its dependence on the point where the 
tensor fields g ab and K ab are evaluated, which depends on F itself. This means that g ab and K ab play the same role in 
the apparent horizon equation as the internal metric of a nonlinear er-model does in its equation of motion. Because 
the coefficients of the elliptic equation contain g ab and K ab as free functions, it appears unlikely that one can prove 
existence of solutions for sufficiently general g ab and K ab . 

C. Geometric characterization of the L 2 operator and spherical harmonics 

In this subsection we introduce in geometric terms some tools that we need later on to discuss spectral methods 
on the 2-sphere. The key idea of any pseudo-spectral method for solving a nonlinear elliptic problem is to subtract 
from the nonlinear one a simple linear elliptic operator that can be inverted explicitly by spectral methods. In our 
problem, the principal part of the operator H acting on h is the Laplacian with respect to the 2-dimensional metric 
m ab induced on the surface F = by the metric g ab - As F — is topologically a 2-sphere, a natural candidate for 
subtraction is the Laplacian L 2 on the round 2-sphere. It can be inverted using the spherical harmonics. 

We could define L 2 on E by first introducing spherical coordinates (r, 9, ip), and then defining its action on a scalar 
/ as the usual combination of partial derivatives with respect to those coordinates, that is, 
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L J f = f,eo + cot Of,e + sin" 2 (11) 

Setting up these coordinates also has the effect of lifting the spherical harmonics from the 2-sphere to all of E, by 
smoothly identifying points on different spheres r = const. 

The minimal geometric structure which allows us to make the same definitions without reference to preferred 
coordinates is a flat background metric g a b on E (independent of the physical metric g a b), together with a preferred 
point O. Let the covariant derivative associated with g a b be V a , and let g ab be the inverse of g a b- We foliate E into 
level surfaces of the scalar field r, where r(p) is the geodesic distance with respect to g a b between the points p and 
O. The vector r a = g ab Vf,r is the unit normal with respect to g a b on the surfaces of constant r. The flat metric g a b 
then induces the metric g a b — V a rV;,r on the surfaces of constant r. This induced metric has a constant curvature 
of r~ 2 , so that r~ 2 (g a {, — V a rV;,r) is a metric of unit curvature on the 2-spheres r — const. We now define L 2 as the 
Laplacian of this 2-dimensional metric: 

L 2 = r 2 (g ab - r a r b )V a V b - 2rr a V Q . (12) 

By direct substitution one verifies that, if g a b is given as 

ds 2 = dr 2 + r 2 {d0 2 + sin 2 9dp 2 ), (13) 



this reduces to Eq. (|il|). Our definition (|12|), however, is covariant, and can be used to define the action of L 2 on 
arbitrary tensors, and in arbitrary coordinate systems. 

For our purposes we characterize the spherical harmonics Y/ m as a set of scalar functions on E with two properties: 
They are orthonormal in the sense that 

mm' 1 

(14) 

Js 

where S is any smooth surface that is star-shaped around r = 0, and where dfl is the measure induced on S by 
r~ 2 {gab — V Q rV;,r). (In spherical coordinates this reduces to the standard measure dfl = smOdOdtp.) From this it 
follows that 

r a V a Y lm =0. (15) 

[In spherical coordinates Yi m — Yi m (0, ip).] We also require that the Yi m are eigenfunctions of L 2 : 

L 2 Y lm = -1(1 + l)Y lm , (16) 

for I = 0, 1, 2, . . . and m = —I, . . . , I. We do not define m as the eigenvalue of L z (8/ dip in spherical coordinates), but 
only use it as a label on the orthonormal basis. This leaves us free to combine Yj m and Yi,- m of the standard complex 
definition to obtain a real orthonormal basis more convenient for numerical purposes. 



III. ALGORITHMS FOR SOLVING THE APPARENT HORIZON EQUATION 



A. The Nakamura et al. algorithm 

We now use our covariant notation for L 2 and the Y\ m in reviewing the algorithm of Nakamura, Kojima and Oohara 
H (from now on NKO) for finding an apparent horizon. NKO characterize 5* by r = h(9, <p) in spherical coordinates, 
and expand h in spherical harmonics: 

h{0,<p) - £ a lmYlm{G,9)- (17) 

1=0 m=-l 

(A finite value of Z max is required in any numerical implementation.) We begin our description of the algorithm with 
the trivial observation that H = is equivalent to 

pH + L 2 h = L 2 h, (18) 

where p is any strictly positive function. In the NKO algorithm, the weight function p is specified by demanding that 
the coefficient of the partial derivative h t gg cancels in the combination pH + L 2 h. (The notation p is ours, not that 
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of NKO. We introduce it here because we want to consider other choices of p later on.) Integrating over the Yi m and 
using ( |l4|) and (|l6|), we obtain 

/ Y l * m (pH + L 2 h)dn = -l(l + l)ai m . (19) 
Js 

NKO now use this equation in an iteration procedure, {ai m }( n > — > {a; m }^™ +1 \ where (n) labels iteration steps, of the 
form 

= - j^i) / s Y *m{pH + L 2 h)^cK}, (20) 

where the right-hand side is evaluated from the {aim} ■ As this formula does not cover aoo, aoo is determined at 
each iteration step by solving 

/ (pH + L 2 h)dn = (21) 
Js 

for ooo (note Yq = const.) 

We now try to understand what makes the NKO method work. For this purpose we express H(h) in terms of the 
flat background derivation V a : 

V a V fc F - ig Crf (V a5cfc + Vf, 5ac - Vcffafe)VrfF 



H= (g ab -s a s b )( \VF\ 



-K ab }, (22) 



where F = r — h(9, ip), and s a is defined by Eq. (g). Because r a V a h — by definition, we have 

L 2 h = r 2 {g ab -r a r b )V a V b h. (23) 

Putting Eqs. ( p2| ) and ( p3| ) together, and keeping in mind that F = r — h, we obtain 

pi? + L 2 /i = M ab V a V b h + W, (24) 

where M ab and W depend explicitly on first derivatives of h, and implicitly on h through the point in S where g ab 
and K ab are evaluated. We have quietly assumed that p does not depend on second or higher derivatives of h, so that 
pH + L 2 h, like H itself, acts on ft as a quasi-linear second-order differential operator. This is indeed the case for the 
p of NKO and the other choices we explore later on. The principal symbol M ab is 

M ab = -p\VF\~ 1 (g ab - s a s b ) + r 2 (g ab - r a r b ). (25) 

The principal symbol of a quasi- linear differential operator does not depend on the choice of derivation, here V a . We 
can verify this for the case at hand. 

We see that M ab is a difference between two projectors: the first one onto surfaces of constant F, and with respect 
to the physical metric g ab , the other onto surfaces of constant coordinate r, and with the respect to the background 
metric g ab . In the trivial case where g ab is conformally related to g ab (conformally flat), and where surfaces of constant 
F coincide with surfaces of constant r (spherical symmetry), one can choose p so as to make the entire tensor M ab 
vanish. In general, one can impose only one condition on its six components. The choice of NKO is, in our notation, 

M ee = 0. (26) 

We prefer a coordinate-independent choice, and impose 

M ab (g ab - V a rV b r) = 0. (27) 

The motivation of either choice is to cancel, as far as possible, that part of M ab V Q Vf, which looks like L 2 . Our choice 
does not introduce a preferred direction within the tangent space of 5, which may be an esthetic more than a practical 
advantage. Solving our condition for p, we obtain 

p = 2r 2 |VF| [(g ab - s a s b )(g ab - V a rV fc r)] ^ = \VF\a, (28) 
where the second equation defines a. 
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Now we recognize an important ingredient of the NKO algorithm, its smoothing property. Putting the individual 
components ai m back together again, we can write ( ^o|) as 

h {n+i) = (tfyi^H + £ 2 )/j(™). (29) 

(This is only formal because of the special role of aoo: L 2 does not have an inverse.) Any iterative algorithm for solving 
an elliptic problem runs the danger of being unstable to the growth of high-frequency numerical noise. Whereas H 
acts on h as a second-order differential operator, thus increasing unsmoothness, pH + L 2 has the L 2 part taken out, 
and therefore creates less high-frequency noise. Moreover, (L 2 ) -1 acts as a smoothing operator. One therefore expects 
to be smoother than h^. This is a necessary property for any iterative algorithm that can converge from a 
rough initial guess without blowing up through high-frequency noise on the way. 



B. Jacobi's method, and stability 

In order to see how the NKO algorithm is related to other algorithms, we rewrite ( p9| ) once more, as 

h (n+i) _ h (n) = (L 2 )-\pH)( n \ (30) 
It is now tempting to go from the discrete algorithm to a continuous flow in an unphysical "time" parameter A: 

dh(9,<p,\) , r2W 



dX 

The NKO algorithm proper, namely 



(LT L pH(h). (31) 



(n+l) (n) _ 1 / TT\(n) 



a lm a lm ~ 1(1+1) KH±± ,lm ' 

is formally recovered from this differential equation by forward-differencing it with respect to A, with a step size 
AA = 1, and inverting L 2 by the pseudo-spectral method (we again disregard the special role of aoo). Other differencing 
methods, such as centered differencing, and using a different "time" step, give rise to obvious alternative algorithms. 
Some of these have been examined by Kemball and Bishop || . Kemball and Bishop also consider different methods 
of enforcing the constraint ( pl| ) on aoo and of coupling it to the iteration method for the other ai m . 

Any flow method can be considered as an example of Jacobi's method. This is the recipe of solving an elliptic 
equation E(h) — by transforming it into a parabolic equation dh/dX = E(h). If E is the Laplace operator, then 
the resulting equation is the heat equation, and Jacobi's method is known to converge. As H acting on h resembles 
— L 2 , one might try the flow 

§ = -H{h). (33) 

We have implemented this numerically in the pseudo-spectral framework and find empirically that its high frequency 
noise blows up unless one chooses a very small step size. 

The origin of this instability is clear from the analogy with the heat equation. The heat equation on S 2 is 

f.x = L 2 f. (34) 

We decompose / into spherical harmonics, as f(9, ip,X) fim(X)Yi m (9, <p). For the spectral components we obtain 
dfim/dX = —1(1 + l)fim- All spectral components decrease exponentially. Discretizing this equation in time, however, 
for example by forward differencing, we obtain//™ +1 ^ = ff^ ' —A\l(l+l)f^ '. This is stable only for |1 — AAZ(Z+1)| < 1, 
that is for AA < 2/1(1 + 1). In all explicit methods, the time step is limited by order of magnitude to 

AA ~ 1 (i~ +iT (35) 

The same limit arises if one discretizes L 2 by finite differencing, where it takes the form AA < (A9) 2 ~ (Aip) 2 . 
Similar stability limits exist for all parabolic equations. The NKO algorithm does not have this instability problem. 
It replaces H = by (L 2 )^ 1 (pH) = as the elliptic problem to be solved, and clearly (L 2 )^ 1 acts a smoothing 
operator that keeps high frequency noise down. An appropriate choice of p makes this even more effective by making 
pH as similar to — L 2 as possible. 



(PH)Z!, l>0, (32) 
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C. Mean curvature flow 



From considering an iterative approach as the discretization of a flow on the space of surfaces, one is led to the 
generalized mean curvature flow algorithm of Tod H and other geometrically motivated flows. Tod proposes deforming 
a trial surface S embedded in E by means of the flow 



d_ 

dX 



= -s a H, (36) 



where s a is again the outward-pointing normal to S. [Tod uses the notation dx l /dX for the left-hand side, but we 
wanted to stress here that (d/d\) a is a vector field and independent of coordinates.] For time-symmetric Cauchy 
data, K a b = 0, we have H = V a s a , which is simply the trace of the extrinsic curvature of S induced by its embedding 
in S, also called the mean curvature. H = is then equivalent to S having extremal area, and s a H is the gradient of 
the area. In this case, mean curvature flow is guaranteed to converge to a surface of H — 0, or extremal area, also 
called a minimal surface. There is an extended literature on mean curvature flow and minimal surfaces jig] . Tod's 
idea is to generalize this method from H = V a s a to H — V a s a — K + K ab s a s b . For K ab ^ 0, this flow is no longer 
guaranteed to converge, but one may hope that it does, as the additional terms are of lower order. 

One essential strength of generalized mean curvature flow is that it cannot move a test surface through an AH, even 
for K ab 7^ 0. The argument is simple fTi| : Assume that the test surface is about to move through the true AH, that 
is, it touches it at one point. At that point both surfaces see the same g ab , K ab and s a . Of the quantities which go 
into the expression (||) for H, only the V a s a differ on the two surfaces. Keeping track of the signs, one sees that the 
test surface must then always back away from the true AH at that point. Therefore, a smooth test surface can never 
cross an AH (although it can approach it asymptotically). This is true not only for generalized mean curvature flow, 
but for all flows of the form (d/dX) a = —s a pH, as long as p is strictly positive. This property allows us to start the 
algorithm on a large surface far out and evolve it inwards, thus making sure we find the true AH. 

We note that Eq. (|36|) does not only specify the deformation of S as an abstract surface, but also identifies any 
point on S with a point on its deformation. That information is not essential to the method, and we get rid of it if 
we define S as the level set F = 0. Then s a is again given by (||). Consider a family of moving surfaces S(X) given 
by F(x z , X) — 0. On the surface F — const, we must have 

dF dF ( d\ a 

dx^^x + {dx) v ^ = °< (3? ) 

and therefore 



dF ( d \ a 

~dx = ~\dxj v " F = saHV * F = \ VF \ H - ( 38 ) 

We note that (^) is a geometric prescription: It specifies a vector field on 5 only in terms of the geometry of S and 
E and the tensor field K a b, independently of how S is parameterized. As wc have just shown, the parabolic equation 
( |38| ) is equivalent to it. We conclude that a flow parameterized by, for example, dF/dX = H, without the factor |VF|, 
does not have such a geometric interpretation, but must depend on F in a more general way than only through the 
shape of its level set S. On the other hand, as H is a scalar function of g a b, s a and K a b (evaluated on S), we can 
replace it by any other scalar and still obtain a flow with geometric meaning. Any flow of the form 

dF 

— = VF x any scalar (K ab , g ab , s a ) (39) 
oX 

is therefore geometric in nature. Such a general equation, replacing H by any function of the curvature of 5 1 , has 
already been given by Osher and Sethian [ |l8[ . 

If we now restrict F to the form F(x l , X) = r — h(6, ip, A), we have dF/dX — —dh/dX. Therefore, any flow of the 
form 

dh 

— = -|V(r- h)\ x any scalar(if afc , g ab , s a ) (40) 
is again geometric in nature. The naive Jacobi method, Eq. (|33|), however, is not. 
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D. "Fast flow" methods 



Before we propose our own AH rinding algorithm, we summarize the strengths and weaknesses of the existing ones. 
We have not discussed algorithms which attack the elliptic problem directly via Newton's method or a minimization 
iteration. Their main drawback, however, is a small range of convergence, that is, they require a very good initial 
guess. NKO is a lot more robust, but the need to treat aoo separately is an important disadvantage. Eq. is by 
no means trivial: Solving it by any iterative method like Newton's method is as computationally expensive as many 
steps of the main iteration loop. Furthermore, for given data (g a b,K a i,) and given a; m , I > 0, there may be several 
roots aoo of this equation, or none. Kemball and Bishop propose investigating each root of the equation separately, 
or if there are none, each minimum of the right-hand side. Clearly this makes the algorithm even more expensive. 

Most importantly, Newton's method for solving ( f2l|) tends to go off into the wrong direction. As an example for this 
problem JlS|] , we consider time-symmetric Cauchy data for the Kruskal spacetime of mass m in isotropic coordinates. 
[This is the special case = of Eq. ([H]) below.] As a trial surface we take a sphere of coordinate radius f centered 
on the black hole. The expansion of outgoing null rays is 

(2r + m) A 

From a mathematical point of view, this example is degenerate in the sense that H = is a reduced from a differential 
to a purely algebraic equation by the spherical symmetry of the trial surface. (There is only one spatial direction, and 
this is the degenerate direction of the elliptic operator.) Nevertheless, any problems that arise in this toy equation 
also arise in a more realistic situation. From a plot of H(f), Fig. 1, we see that for f > 1.87m Newton's method 
wanders off to infinity, and for r < 0.13m goes towards f — 0, instead of finding the zero at f = m/2. All algorithms 
which use Newton's method, or a minimization method using derivatives, for any or all of the ai m , that is, the direct 
elliptic algorithms, share this problem. 

The curvature flow method is sensitive only to the sign of H , not its derivative. Applied to this problem, it goes 
towards smaller r for positive H, and towards larger r for negative H, and always finds the apparent horizon. We 
have already seen that it cannot accidentally walk through the AH. In these two properties lies its robustness. Any 
flow with pH instead of H on the right-hand side, where p is strictly positive and a scalar, shares the fundamental 
advantages of the generalized mean curvature method: a trial surface far outside the apparent horizon always moves 
in, and can never accidentally cross the apparent horizon. As we have seen, however, flow methods are slow because 
explicit discretizations in time of parabolic methods require very small time steps for stability. An implicit time 
step may be possible, but introduces a new elliptic problem which a priori is not simpler than the underlying elliptic 
problem one wants to solve. 

We are not bound to geometrically motivated flows, however. Instead, we heuristically consider all flow methods as 
variants of the Jacobi method for solving H — 0. Then we are free to combine the best features of the curvature flow 
and NKO methods. From curvature flow we would like to keep the properties that aoo is not treated specially, and 
that the change of aoo should be proportional to -ffoo? not di/oo/^oo- From NKO we would like to adopt the idea of 
subtracting and then inverting L 2 , in order to suppress high-frequency noise. This leads us to the following family of 
flows: 

^■ = -A(l-BL 2 y 1 pH, (42) 

where A and B are free positive constants, and where p is a strictly positive weight depending on h through at most 
first derivatives. The differential operator 1 — BL 2 is invertible, with positive eigenvalues, for B > 0, and for B > 
its inverse is a smoothing operator. When we discretize in A, we can absorb AA into A. For simplicity we also restrict 
ourselves to forward differencing. In spectral components, we obtain 

a <™ _ a *™ " ~1 + BZ(I + 1) {pH)lm ■ (43) 

For p we consider three choices: p = 1 ("H flow"), p = \WF\ ("C flow"), and p = \VF\a with a defined in @ ("N 
flow"). With A > and B = 0, H flow is the Jacobi method, and C flow is the curvature flow method. N flow 



formally becomes the NKO method [compare Eq. (32)] in the limit A — B — > oo. The limit is singular because the 
NKO method is not a flow and has to update the component aoo separately. 

For determining the optimal values of A and £?, it is convenient to reparameterize them with new parameters a 
and (3 as 
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A and B now scale with l max in such a way that we expect the optimal values of a and fj to be independent of the 
value of Z max . a parameterizes an ^-independent contribution to the effective step size of a/[Z m ax(imax + 1)], while j3 
adds an ^-dependent speedup which is zero at I — i max and increases to [3 at I — 0. 

It is clear that the fast flow methods have the potential to be much faster than curvature flow, while still being 
numerically stable and robust against bad initial guesses. They are not really flows of the form ( f4(i| ) because they are 
not local. In some situations, the effective weight p can become negative on parts of the surface, and in these situations, 
the "fast flow" can move through the true AH. Fast flow methods should be considered as (good) compromises between 
the robustness of curvature flow and the speed of NKO. Furthermore, one can trade robustness for speed by increasing 
/3, and vice versa, and so adapt the algorithm to the situation. We have obtained the best results using N flow, with 
a = 1.0 and (3 = 0.5 in (|44|). We note that a = 1.0, for any /3, means that the algorithm treats high frequency 
components like the NKO algorithm does. 



IV. NUMERICAL IMPLEMENTATION OF PSEUDO-SPECTRAL APPARENT HORIZON FINDERS 



The algorithm we suggest in this paper is formally defined by Eq. (|43|) with p defined by (28) and A ~ B ~ 0.5. 

In order to implement this or any other pseudo-spectral algorithm, we need to calculate the spectral components 
(pH)i rn from the spectral coefficients ai m . In this section we give details of an algorithm for doing this, given g %3 , Kij 
and gij t k on a Cartesian grid. We expect that there is scope for increasing the speed and reducing the discretization 
error in this low- level part of the algorithm, without changing the top-level part given by (0) . 



A. The background structure 



The parameterization of the surface S through spherical harmonics and the introduction of the differential operator 
L 2 require the introduction of a flat metric g a b- We do this by introducing auxiliary Cartesian coordinates x l — /'(x- 5 ), 
and then setting the components of g a b in the coordinates x 1 to be 8tj. The corresponding metric derivation V a is 
then d/dx l , and r 2 = 8ijX l xf In these coordinates L 2 is given by the expression 

While more complicated choices are possible, we define 

x l EE x l - xl, (46) 

where x^ are the Cartesian coordinates in which the Cauchy data are presented to our algorithm. The freedom to 
shift the origin r = around is necessary because any trial surface will have to be star-shaped around r = 0, that is 
around Therefore we have to make sure that x l — xl, is inside the AH. 



B. Calculating the Yj m 

We need to calculate the Yi m (x l ) and their first two partial derivatives for arbitrary (x l ). Speed is important, 
because our algorithm spends most of its time in these calculations. The standard spherical harmonics are 

y m = p™( cos e)e lmv , (47) 

where the P/ 71 are associated Legendre functions times a constant depending on I and m. Instead of the complex Yi m , 
we introduce the real basis Yi m as 

?j° = P i °(cos0), F; 1 '" 1 = V2P, |m| (cos(9)cos7n(A = %/2P, |m| (cos 6») sinm^. (48) 

The real Y\ m obey the same conditions ( |l4| - |l6|) as the standard complex Yi m , but they are not eigenfunctions of 
L z = d/d(p. At each point x 1 = [x, y, z) we calculate 

cos9 = z/r, sin.0 = p/r, cosip = x/p 1 sm.cp = y/p, (49) 
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where 



r = \J x 2 + y 2 + z 2 , p = \J 1 x 2 + y 2 . (50) 

No explicit evaluation of trigonometric functions is required. Then cos mtp and sin my are calculated as polynomials 
in cos ip and sin ip from the recursion relations 

cos rrap= cos(m — 1) ip cos ip — sin(m — l)(p sin ip, 

sinnup— cos(m — l)ip snap + sin(m — l)ip cosy. (51) 
The P\ m are given explicitly for m — I by 



H(axe)= ^^m ( _ sta9) «, (52) 



and for < m < I — 1 are calculated from the recursion relations 



pr (cost 21 1 



V2l^l oosfli^cosfl) - v 2 ^_ 3 P ; m 2 (cos( 



(l-l) 2 ~m 2 



(53) 



(They are not needed for m < 0.) In order to calculate the first and second partial derivatives with respect to x, y 
and z, we calculate the partial derivatives of Y\ m with respect to 9 and ip, and those of 9 and ip with respect to x, y 
and z, and explicitly code all terms arising from the chain rule. The derivatives of P; m (cos0) with respect to 9 are 
obtained recursively after differentiating Eqns. (^2j) and (|5^). The relations (|l5|) and ([l6]) are then obeyed to machine 
precision by the numerically calculated quantities. 

We are aware of two other algorithms for calculating Yi m (x l ) and their first and second partial derivatives. The 
algorithm of Baumgarte et al. calculates them recursively as polynomials of the r^. We have coded this algorithm 
directly from the detailed formulae in [ pT[ , and find that it scales in time as Z^ax an d i n storage requirement as Z^ax- 
The NCSA/WashU apparent horizon finder Jl2| , p^ | does not calculate the Yj m , but a related basis of smooth functions. 
This basis is not orthogonal, and it is not independent of r at constant 9 and ip. For the NCSA/WashU algorithm these 
properties of the basis functions are not necessary. The calculation of this basis scales as approximately Z^ax in time, 
and as Z^ax m storage p0| . In common with both algorithms, ours is recursive, and does not require trigonometric 
function evaluations. The difference is that it breaks up the Yi m into the product of a function of 9 times a function 
of ip. In consequence it scales Z^ ax in time (it is faster already for Z max = 2), and as /^ ax in storage requirement. This 
optimal scaling comes at a price: the algorithm breaks down on the axis x = y = 0, where cancellations between the 
9 and ip dependent factors in the analytic expressions fail to take place numerically. In practice, one can evade the 
problem by moving any collocation points that come very close to the z axis a small distance away from it, resulting 
in a small error at that point, and a negligible one in the integrals over S. Incorporating the cancellations into the 
code properly requires mixing the 9 and tp dependency by going through an intermediate, over-complete basis called 
"symmetric trace-free tensors" , which is precisely the approach of Baumgarte et al. 



C. Interpolating and integrating over S 

We need to discretize the integral f„ dQ. We take as collocation points on S all those points where 5* intersects a link 
of the three-dimensional Cartesian grid. A link is the straight line between two neighboring points on the Cartesian 
grid. The links which intersect 5* are those on which F changes sign. The advantage of this choice of collocation 
points is that one only needs to interpolate in one dimension. Furthermore the number of collocation points on S 
scales with the number of nearby points on the Cartesian grid, that is, with the available numerical information. 

To find the surface F = the algorithm calculates F on all Cartesian grid points. For this it needs the Yi m on all 
grid points. Although these are required again and again, present technology does not allow us to store ^ max (^ m ax + 1) 
3D arrays for reasonable Z max , so that they have to be recomputed each time. Then the algorithm flags all links on 
which F changes sign. Both operations scale as TV 3 , where N is the linear grid size. We determine by inverse linear 
interpolation where on the link the intersection point is, then interpolate g lJ , gij,k and i^y to the intersection point 
by cubic interpolation. We calculate Fj and F^ directly at the x 1 of the intersection point. For this purpose we need 
the r i and the Yi m ^ and Yi m j. 

The integral J s dfl is now approximated by the sum 
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fdn-iTT^-, (54) 

where the sum is over all collocation points. Let s l be the unit normal on S with respect to the flat background 
metric g a b, 

s l = [6 kl (r, k - h ik )(r, t - h,0] ~ 1/2 <F(rj - h tj ). (55) 

The integration weight w is then given by 

_ _ dtt dA dtt_ r^£_ dN_ _ js 1 ! |s 2 | |s 3 | 

W ~dN~dAdN' dA~~' IX ~ Ax 2 Ax 3 + Ax 1 Ax 3 + Ax 1 Ax 2 ' ( ' 

where Ax 1 are the grid spacings on the three axes. Here d£l is the solid angle with respect to the flat metric g a b around 
r = 0, that is d£l = shi9d9d<p, dA is the surface element on S induced by g a b, and dN is the number of intersections 
of dA with grid links. Note that the expression for dN /dA models the anisotropy of the Cartesian grid in an explicit 



sum over the three grid directions. The sum (54) is a good approximation to the integral when s l changes little from 
one collocation point to its neighbors. 

As a test of the discrete approximation to integrals over the surface, we calculate the overlap integrals (]l4|) numer- 
ically. Let us combine the indices I and m of the spherical harmonics into one index n. The numerical approximation 
to the symmetric matrix A nn i = J Y n Y^ dil is not exactly equal to the unit matrix, because of the finite number of 
collocation points. 

One could arrange the weights w such that A nn i comes out right for a given set of collocation points, but that would 
require putting the collocation points in fixed, special positions with respect to 9 and ip, for example on a square grid 
in 9 and ip. In our algorithm, however, we let the position of the collocation points be dictated by the underlying 
Cartesian grid x l , and rely on a number of collocation points much larger than the number rt max = ^maxGmax + 1) of 
basis functions in order to keep the error down. Table | shows how the error in A nn > increases with I for a surface 
S created under realistic conditions by the apparent horizon finder. In practice, the size of the 3-dimensional grids 
is limited by the available computer storage, so that we have to choose Z max small enough for the spectral error to 
remain small. 

We can reduce this error in the following way. Let us denote by H„ = H n {a n ') the true spectral components of 
the function H on the surface parameterized by the expansion coefficients a n . Let H n be their numerical, slightly 
incorrect value. Clearly we have 

oo 

H n = ^ A nn >H n i. (57) 

n' = l 

Without much additional numerical work, we can calculate a finite square piece of the infinite matrix A when we 
calculate H n up to n max . Let B be the inverse of that finite part of A. Then we have (for n < n max ) 

n mSLX n max / oo \ oo 

H„ = ^ B nn iH n i = H n + ^ B nn i I ^ A n > n "H n 'i ) ~ H n + ^ A nn >H n i. (58) 

n' — l n' — l \n' f -n max + l / n' — n max + l 

In H n the unwanted aliasing among the low (n < n max ) frequencies has been eliminated, and the remaining deviation 
from the true value H n comes only from the aliasing of high frequencies to low ones. One would assume this to be 
a better approximation to H n than the H n in normal situations. In practice, however, this assumption is difficult to 
test, as no cheap estimate of the error H n — H n is available. 

Still, there are some indications of the remaining error in the H n : The spectral components of r, the r„, are by 
definition identical to a n . We find that the f„ are much closer to a n than the f„, but do not converge to them. 
The remaining error can only be due to the fact that the collocation points do not lie exactly on the true surface 
S parameterized by the a n , due to the interpolation used to find them. This failure to find the true surface S is 
the only source of error for the r„, but appears to be also the dominant source of error for the H n , or any other 
nontrivial function on S. In practice we proceed as follows: We use H n as our best approximation to H n . We monitor 
convergence of the final result a n of the AH finder with Z max and the grid spacing of the underlying Cartesian grid. 
We also monitor |f„ — a n \. Finally, and perhaps most importantly, we find that the algorithm using H n converges 
better, and its error is considerably reduced when tested against data for which the apparent horizon is known in 
closed form. Therefore we always use the H n and other hatted quantities in the algorithm. 

In order to estimate the quantity 
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/ H 2 dn = Y,H 2 n , (59) 
J s n =i 

which indicates to what precision the algorithm has found the apparent horizon, we use the two numerically available 
quantities 

(ff rms ) 2 = = Y,Y, A nmH n H m , and \H\ 2 = ]T H 2 n . (60) 

n—l m— 1 n— 1 

After this work was carried out, we became aware of a different, perhaps more efficient algorithm for going back 
between a function of 8 and </> and its spherical harmonic components pl[ . One puts a grid on S which is rectangular 
and equally spaced in 9 and (p, and then uses fast Fourier transforms in 9 and (f>. In a second step, one has to discard 
those linear Fourier components which are not sufficiently regular at the poles 9 = 0, tt, which is rather complicated. 
In order to evaluate gij etc. at the collocation points required now, one has to interpolate in three dimensions from 
the given Cartesian grid, instead of in one dimension. Nevertheless, there may be scope for a more efficient algorithm 
here. 



V. TESTS 



A. Brill-Lindquist data 

The NCSA/WashU algorithm appears to be the only one to have been tested on numerically evolved data Jl3| . 
Tests that use data given in closed form avoid interpolation from numerical data on a grid, which poses an additional 
source of numerical error and even instability in realistic situations. In the present paper we test our algorithm with 
data given in closed form, but passed to the algorithm only on the points of a numerical grid of realistic size. The 
input into the code are the numerical values of the inverse metric and extrinsic curvature components and of the first 
partial derivatives of the metric components on the grid. A performance test with data derived both from numerical 
initial data algorithms and numerical time-evolutions is left to a future publication. 

As a first test of the complete AH finder, we use Brill-Lindquist time-symmetric initial data. For two black holes, 
these are 

'H^*^^)^ <61) 

The generalization to N black holes is clear. 

We begin with a single black hole, where the AH is known explicitly: it is a coordinate sphere of radius m/2. There 
are a number of possible convergence criteria for the iterative algorithm, none of which fits all possible situations. 
One such criterium is fl rms > 2\H\. [These measures were defined in (|60|).1 This means that the residual of H(9,(p) 
is mainly in the high frequencies that we do not resolve. Table [□] shows the performance of the algorithm for this 
convergence criterium. We chose a grid spacing such that there are roughly 16 grid points across the interior of the 
AH. By the standards of a 3D single grid numerical relativity code on a current supercomputer that is already as 
much resolution as one can hope for. We chose / max = 6, which is roughly the optimal value for that resolution. The 
initial data are always aoo = 0.8, while the horizon radius is 0.5. (For convenience, we use ai m rescaled by a factor of 
a/47t, so that aoo is the average coordinate radius of the AH.) 

We have varied the offset of the center of the spherical harmonics from the center of the AH. Ar is the error in 
locating the AH in coordinate space. It is calculated directly in those points where the AH is collocated by the 
algorithm, not from the <z/ TO . The result is roughly independent of the direction of the offset. We see that if the 
surface is very eccentric around the origin of coordinates, precision suffers. Fortunately, there is a simple remedy: If 
the dipole moments aj^ 1 are large, the algorithm automatically uses them to obtain a better value for the origin Xq 
of coordinates, and restarts. This source of error is totally eliminated by the procedure. The same remedy applies if 
the surface touches the origin of coordinates at any stage during the flow. 

The next test is Brill-Lindquist data for two uncharged black holes of equal mass. We can position the two centers 
x\ and x-i so that the metric is symmetric with respect to the x = 0, y = and z = planes. This allows us to work 
numerically on an octant of the full grid, and save time and storage. The situation is in fact axisymmetric, but the 



code does not know that. In Fig. VI all points on the discretely represented surface are plotted, giving coordinates z 
versus p = y/ 'x 2 + y 2 . The fact that they all fall on one curve shows that the code represents an axisymmetric surface 
well in spite of the underlying Cartesian grid. 
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In the data (|6l|) one can always find two minimal surfaces surrounding x\ and 12- If x\ and X2 are close enough 
together, there is a third minimal surface surrounding both of them. Determining the maximal separation at which 
this happens is not an easy test. Assume that the two centers are just far enough apart that there no longer is a 
common horizon. By continuity there will still be a smooth surface on which H is small, but not zero, everywhere. 
Numerically, this cannot be distinguished from a true horizon. 

In the test, the two black holes have equal mass parameters mi = 7712 = 1. The total ADM mass is 2. We look for 
both inner and outer surfaces. In Table III we show, for the same numerical parameters, the root-mean-squared value 
of H on the trial surface after our algorithm has stopped, against the (coordinate) separation d of the two centers 
x\ and X2 ■ The axisymmetric numerical algorithm of Brill and Lindquist does not find an outer minimal surface for 
d > 1.56. From our calculations we can say with confidence that the limit lies between 1.4 and 1.6. We should stress 
again that this precision is limited by the resolution of the Cartesian grid on which we give the Cauchy data. An 
algorithm specialized to axisymmetry could of course determine this limit with much higher precision. 



B. Kerr data in Cartesian coordinates 

In order to test our code on analytic data with nonvanishing extrinsic curvature, we consider Kerr data. Cauchy 
data for the Kerr spacetime have been given by Brandt and Seidel . We transform these to Cartesian coordinates 
by defining x t x k Sij — f 2 , where f is the radial coordinate which generalizes the isotropic radial coordinate f for the 
Schwarzschild spacetime. For testing our algorithm it is useful if we do not restrict the angular momentum vector, or 
symmetry axis, to the z-axis, but give its direction as a unit vector ri 1 , n l n 3 8ij — 1. The transformed expressions are 



gij = A(6ij + BviVj), g i3 = A 1 f S ij - l - -^—^ v*^ j , =v i w i +w i v i) (62) 

v*= e ijk njXk, w i = Cx l + D(n l - cos9x i ), (63) 

where cos# = n k Xk/r, where all indices are moved with <5y, and where the coefficients are 

A=p 2 f- 2 , B = (p 2 + 2mr)a 2 p~ 2 f~ 4 , (64) 

C= [(r 2 + a 2 ) 2 - Aa 2 sin 2 9] ~ 1/2 am [2r 2 (r 2 + a 2 ) + p 2 (r 2 - a 2 )] p" 3 ?" 4 , (65) 

D= [{r 2 + a 2 ) 2 - Aa 2 sin 2 9] ~ 1/2 2a 3 mrA 1 / 2 cos^f" 4 , (66) 

2 2 

oono„.n o 7TI 0/ , > 

p = r + a cos V, A = r — 2mr + a , r = m + r H = — . (67) 



The apparent horizon is the coordinate sphere f = \]m 2 — a 2 /2. For a — these data reduce to Brill-Lindquist data 
for a single black hole. While K a f, does not vanish, the data are still special in that the two contributions to H, 
V a s a and m ab K a b, vanish separately on the AH. We have tested our algorithm on these data for different ratios of 
a/m, different offsets x l between the center of the black hole and the center of spherical harmonics, and for different 
orientations n l of the black hole symmetry axis relative to the spherical harmonics. The results are essentially the 
same as for the single Brill-Lindquist black hole, giving an indication that the presence of the extrinsic curvature term 
does not make a qualitative change to the performance of the algorithm. 



VI. CONCLUSIONS 



Numerical general relativity requires a fast and robust algorithm for finding apparent horizons in Cauchy data 
without symmetries in three dimensions given on a grid. In this paper we have described a new apparent horizon 
finder algorithm which appears to be as fast but more robust than its best predecessor. 

We began from a general classification of possible approaches to the problem. Any approach which poses the 
problem as a nonlinear elliptic equation on a topological two-sphere and then attacks that equation directly will fail 
unless provided with a very good initial guess, because the problem is nonlocal in nature. While we have disregarded 
such approaches here, they will be ideal as a second stage whenever the apparent horizon needs to be determined to 
high precision. We concluded that for robustness the best algorithm is probably the (generalized) mean curvature 
flow suggested by Tod 0, where an arbitrary initial surface evolves in an unphysical "time" towards the apparent 
horizon, turning the problem from an elliptic into a parabolic one. The algorithm is guaranteed to converge at least 
for time-symmetric (K a b = 0) data, and we have argued that it must be at least very robust also for K a b =/= 
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data. Unfortunately, numerical implementations of this algorithm face a numerical stability problem common to all 
parabolic equations, which make them slow, and increasingly so with increasing resolution, in practice. 

This stability or speed problem is not present in the algorithm of Nakamura et al. (NKO) B. It is motivated by a 
standard way of solving nonlinear elliptic problems numerically, namely subtracting a simple linear elliptic operator 
from the nonlinear one, inverting it by pseudo-spectral methods and iterating. Here we have thrown more light on how 
NKO works, by making explicit the background metric it introduces, and by characterizing the iteration procedure as 
a specific finite differencing, in unphysical time, of a parabolic problem. This parabolic problem itself is the singular 
limit of a certain family of flows which are governed by a mixture of the physical geometry of the Cauchy data and an 
unphysical background geometry. Tod's flow is a different limiting case of that family, one in which no background 
metric appears. 

Once we have recognized the existence of a continuum of possible algorithms between Tod and NKO, it is plausible 
that an algorithm somewhere in the middle of the continuum may be better than the extremes. By trial and error, 
we have determined the optimal member of the family of algorithms. This intermediate algorithm evolves the high- 
frequency components (the fine details) of the trial AH essentially like the NKO algorithm, but it evolves the low 
frequency components (the rough shape) by a variant of generalized mean curvature flow. We therefore call it "fast 
flow". 

We have given details of a numerical implementation of the pseudo-spectral methods which are needed for imple- 
menting both the original NKO and our new algorithm. Such details have not been published before. It should be 
stressed that the formal analysis of the algorithm in Sect. Ill is independent from its implementation in Sect. IV, 
and there may be different and more efficient implementations. 

We have not made direct performance comparisons with other algorithms, and the tests we have described are 
viability rather than performance tests. Nevertheless, we anticipate the following: 

By construction, the new algorithm is as fast as that of NKO: The iteration steps are very similar, and there is the 
same small number of them. NKO, however, updates aoo (the overall radius of the trial AH) by a special procedure. 
According to how this is done jfj), the additional overhead may be large. More importantly, the separate update of 
aoo has the potential to reduce the robustness of NKO: Eqn. ( |2l|) may have several solutions, in which case all should 
be investigated, or none, in which case minima should be investigated instead ||. This requires some decision-taking, 
which will be hard to automate, or instead an infinitely branching search. We have also argued that zeros of eqn. (|2l| ) 
are hard to find. Either NKO or fast flow should be far more tolerant of initial guesses than the elliptic methods. 

The method of choice for robustness and elegance is clearly Tod's mean curvature flow. The only question here is 
speed. We have argued that as a parabolic method this would be slow in the possible implementations known to us, 
but a quantitative comparison with the implementation of Pasch |j) would be interesting. 

The only existing algorithm directly accessible to the author is that of the NCSA/WashU group |l^,|u|. Direct 
comparisons are planned in future realistic applications. The new algorithm seems to be more robust, though: it 
made the transition from two AHs to a single one in the family of Brill-Lindquist data summarized in Table 1 without 
external input. This ability will be crucial for applying AH boundary conditions in the merger of two black holes when 
the computer code is on its own during a very large run. The NCSA/WashU code has not been tested on a similar 
sequence of analytic initial data, but in some situations involving evolved black holes it presently requires some care 
in finding the correct horizon 23 1. 

Finally, the source code of the new spectral AH finder will be published early in 1998 in conjunction with the 
" Cactus" numerical relativity infrastructure [E4| . 
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TABLE I. Maximal deviation of the overlap matrix A nn i from the unit matrix as 
integration surface S and / max 



function of the linear grid size, the 



Surface parameters 



grid size 



8 



16 



aoo = 1-0 

a 00 = 1.0, ai,_i = 0.4 
a o = 1-0, aio = 0.4 



16 

32 
64 
32 
64 
32 
64 



0.053 
0.016 
0.011 
0.039 
0.047 
0.052 
0.050 



0.159 
0.037 
0.026 
0.069 
0.047 
0.087 
0.056 



0.284 
0.069 
0.040 
0.178 
0.053 
0.274 
0.075 



TABLE II. Root-mean-square residuals of H and maximal errors in the position of the numerically calculated AH, for 
Schwarzschild data offset from the coordinate origin. The coordinate radius of the AH is 0.5. 



Offset 



iterations 



Hrn 



9 x 10" 4 
9 x 10 -4 

1 x 10~ 3 

2 x 10~ 3 
6 x 10~ 2 



( Ar) mtt3; 
7 x 10" 4 

7 x 10~ 4 

8 x 10~ 4 

1 x 10~ 3 

2 x 10" 2 



0.0 
0.1 
0.2 
0.3 
0.4 



10 
10 

11 

12 
10 



1G 



TABLE III. Root-mean-square residuals of H on the inner and outer numerically calculated minimal surface in Brill-Lindquist 
data for two black holes of equal mass. "No convergence" is the unaided return status of the algorithm. It means that the 
residual value of H given in brackets is not due to a lack of numerical resolution. 



Separation inner /f rlns outer ff r]lls 

0.0 only one surface 1.9 x 10~ ; ' 

0.4 8 x 10~ 2 1.8 x 10~ 5 

0.8 9 x 10~ 5 1.8 x 10~ 5 

1.2 3.0 x 10~ 4 1.5 x 10~ 4 

1.4 2.6 x 10~ 4 2.0 x 10" 3 

1.6 2.8 x 10~ 4 no convergence (3.0 x 10 -2 ) 

1.8 2.4 x 10~ 4 no convergence (3.9 x 10 _1 ) 

2.0 7 x 10~ 4 (not attempted) 
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-0.4 



FIG. 1. Plot of the horizon function H(r) versus r, in units of the black hole mass m, as given in Eq. 




-0.2 1 1 1 1 1 1 1 1 1 1 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

FIG. 2. Shape of the AH in the axisymmetric, z-reflection -symme tric situation. The algorithm assumes x, y and z-reflection 
symmetry, but not axisymmetry. The plot shows z versus y/x 2 + y 2 for all grid points on the AH in one octant of the full grid. 
The small half circles are the inner horizons, for a separation of d = 0.4, 0.8, 1.2 and 1.4 of the two black holes, from bottom 
to top. The large quarter circles are the outer horizons, from right to left, or bottom to top. 
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